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Abstract 

We present three dimensional fluid simulation results on the temporal evolution and nonlinear 
saturation of the magnetic curvature driven Rayleigh- Taylor (RT) instability. The model set of 
coupled nonlinear equations evolve the scalar electric field potential cp, plasma density n and the 
parallel component of the magnetic vector potential ip. The simulations have been carried out in 
two limits, (i) a low resistivity case in which RT is the only linearly growing mode, and (ii) a high 
resistivity case where the drift wave is unstable and for which the magnetic curvature parameter is 
set to zero to ensure the absence of the RT growth. Our simulations show nonlinear stabilization in 
both these limits. The stabilization mechanism is similar to that observed in earlier two dimensional 
simulations, namely the generation of zonal shear flows which decorrelate the radially extended 
unstable modes. However the nature of the saturated nonlinear state in the 3d case differs from that 
of 2d in some important ways such as by having significant levels of power in short scales and by the 
presence of electromagnetic fluctuations. Though, in the linear regime the electromagnetic effects 
reduce the growth rates, in the nonlinear regime their presence hinders the process of stabilization 
by inhibiting the process of zonal flow formation. Thus the parameter regime for which nonlinear 
stabilization takes place is considerably reduced in three dimensions. 



I. INTRODUCTION 



The magnetic-curvature-driven Rayleigh- Taylor (MCD-RT) model is a useful paradigm 
for the study of long wave length nonlinear shear flow patterns such as zonal flows and 
streamer structures |l| . These coherent potential patterns are believed to play an important 
role in determining the turbulent transport of matter and heat across field lines in magneti- 
cally confined plasma devices such as tokamaks, stellarators etc. □□□□□□ Q lialiil La- 
in the past few years a large number of investigations have been devoted to the elucidation 
of their characteristics and in understanding their contributions towards transport processes 



In the tokamak context, theoretical studies have mainly 
centered around the drift-wave model or its several variants such as the the ion temperature 
gradient (ITG) mode, the electron temperature gradient (ETG) mode .etc. 
The basic physics underlying the formation of these nonlinear structures is the onset of 
a long scale modulational instability arising from the nonlinear parametric interaction 
of a large number of short scale fluctuations; this mechanism is in essence common to a 
number of earlier model calculations of strong plasma turbulence including that of the 
coupled interaction of Langmuir waves and ion-acoustic excitations llol. In a recent work 

n 

[1| the MCD-RT model was used to explore a very fundamental question pertaining to the 
formation of these nonlinear patterns, namely what determines the symmetry of the final 
state - whether one gets zonal flows or streamer formations Q. The study highlighted the 
role of dissipative processes (namely the strength of the dissipation parameters - the flow 
viscosity /i and the diffusion coefficient D) in influencing the evolutionary path of the system 
towards a particular final symmetry state. Using extensive numerical simulation data it was 
possible to construct a consolidated "phase diagram" in D — \x space which showed that 
low dissipation favoured the formation of zonal flows leading to saturated stationary states 
whereas high dissipation led to formation of radially elongated streamer flow patterns. The 
primary impact of the dissipation was on the initial evolution of the short scale fluctuations 
and the consequent distribution of their spectral power. High power in the short scales 
(due to weak dissipation) at the initial stages led to stronger nonlinear generation of zonal 
flows whereas for high dissipation the streamers gained an upper hand and zonal flows 
were subdominant. We believe that such a unified and consolidated approach could be 
useful in gaining understanding of the evolution of similar nonlinear structures in tokamak 
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transport models as well. Motivated by such considerations and also keeping in mind the 
fact that the RT model in a generic sense can provide valuable insights in a number of other 
experimental scenarios (e.g. currentless toroidal devices, ionospheric spread F irregularities 
etc.) we extend and further develop explorations on the nonlinear dynamics of this model. 
The earlier results reported in ^, Q], were based on two dimensional simulations. In this 
paper we report a major extension of this model by making it fully three dimensional and 
also by including electromagnetic effects through the contributions of magnetic fluctuations. 
In terms of basic physics the extended model now introduces linear and nonlinear coupling 
to shear Alfven modes through finite ku effects so that we now have a set of three coupled 
nonlinear equations that evolve the scalar electric field potential 0, the plasma density n 
and the parallel component of the magnetic vector potential tp. We also have an additional 
dissipation parameter in the form of the resistivity coefficient rj s which makes the drift wave 
branch linearly unstable in certain parameter ranges. We continue to explore the same 
fundamental issues in this generalized model, namely the existence of saturated nonlinear 
states, their characteristics and the factors that influence their formation. Our approach is 
primarily numerical and we present extensive simulation results from our model equations 
to provide answers to the above issues. On the question of the existence of saturated states, 
past results from numerical explorations of the drift wave and allied models have been 
somewhat equivocal and have indicated that electromagnetic effects tend to inhibit zonal 
flow development. In our numerical simulations of the fully electromagnetic RT model 
we find that saturated states still continue to, exist although in a restricted parameter 
domain. Comparison with the earlier two dimensional results show a similarity in the 
saturation mechanism, namely through the excitation of zonal flows. However, there are 
significant differences in some of the characteristics of the saturated states. Thus the 
three dimensional nonlinear states are found to possess a significantly higher power level 
in short scales as compared to their 2d counterparts. Another significant difference is 
that the spatial structures of the potential and density fluctuations do not develop any 
functional correlations. In other words, unlike in the 2d case, the density evolution does 
not slave itself to the potential evolution. Our findings on the effect of electromagnetic 
fluctuations on RT turbulence are similar to those of drift wave models namely that 
they tend to inhibit the formation of zonal flows and thereby to considerably restrict the 
parametric domain of nonlinear stabilization. To highlight the role of the third dimension 
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we have also explored a simplified limit of our model that can be termed as a three 
dimensional electrostatic model. In this simplified model the role of k\\ and the resistivity 
factor rj s are more transparent and we discuss their influence in the formation dynamics 
of the saturated states. In this limit we also demonstrate the stabilizing influence of the 
secondary Kelvin-Helmholtz instability 0] in controlling the unlimited growth of streamers. 

The paper is organized as follows. In the next section we present our generalized model 
equations and and discuss its characteristics including its relation to the previous 2d model 
equations. We also derive a simplified limit of two coupled equations representing the three 
dimensional electrostatic model. Section III is devoted to delineating the properties of the 
linear modes of the model. This is done through numerical and approximate analytic solu- 
tions of the appropriate dispersion relations in various limits. This analysis also highlights 
the role of the various dissipation parameters and the parallel wavelength in the linear evo- 
lution stage of the system. We next present our nonlinear simulation results on the 3d 
electrostatic model in section IV and compare and contrast them with past 2d electrostatic 
results. The full electromagnetic simulation results of the generalized model are presented 
and discussed in section V. The paper ends with a summary of our main results and some 
concluding remarks in section VI. 



II. MODEL EQUATIONS 

The governing equations for the generalized model of the magnetic-curvature-driven 
Rayleigh- Taylor instability are derived along similar lines to that adopted for the previ- 
ously investigated two dimensional set of equations We use the fluid equations of 
continuity and momentum for the electrons and ions along with the quasi-neutrality con- 
dition viz. Vj_ • J± = — V|| J\\. In addition we close the set by the Ohm's law which is the 
parallel component of the electron momentum equation without the electron inertia term. 
We use a slab representation in which the radial coordinate is represented by x, the poloidal 
by y and the toroidal by z. The effect of curvature in equilibrium magnetic field is modeled 
by a x dependent toroidal field B eq = Bq(1 — x/R)z and the radial gradient of the equilib- 
rium plasma density is represented by rio(x) = riooexp(—x/L n ). We take the ions to be cold 
(Tj = 0) but retain a finite electron temperature T e . As a result of this the ion drift (flow) 
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in equilibrium is zero but there is an equilibrium electron diamagnetic drift in the direction 
perpendicular to the magnetic field. Unlike the 2d model we now retain electromagnetic per- 
turbations ( B = B eq + z x Vip), where the perturbed magnetic field fluctuations are assumed 
to arise only from magnetic field line bending perturbations. Thus the vector potential has 
only a z component i.e. A = —ipz. The electric field is given by E = — V0 + (l/c)dip/dtz. 
The scalar and vector potentials are finite only for the perturbations. The total density is 
N = no + n\, thus logN = —x/L n + log(l + ui/uq) = —x/L n + n. The unit vector parallel 
to the magnetic field is given by e» = Bj \ B \= z + (z x Vip)/B . From standard pertur- 
bative expansions, the perpendicular component of the electron momentum equation yields 
the usual transverse flow velocities, namely the E x B drift and the diamagnetic drift. The 
ion equation similarly leads to a E x B drift term as well as a polarization drift contribution 
in the perpendicular direction. These are substituted in the subsequent order equations to 
obtain the nonlinear evolution equations. We express the equations in a dimensionless form 
by normalizing the various physical quantities as follows. The density is normalized by noo, 
the electrostatic potential <ft by T e /e, time by Qi and length by p s = c s /fij. The vector 
potential ip is normalized by Bop s . Our generalized set of model equations then consist of 
the following three coupled equations for the plasma density n, scalar electrostatic potential 
4> and the parallel (to the equilibrium magnetic field B eq ) component of the vector potential 

(J OT1 —* — ♦ I (J —> —* 

— vV + v g — + zxvip ■ vvV - vl { ^viv> + zx • vvl^ = /ivV (2) 

at ay [oz J 

^ + - + y nj~ -ixVf %-n)= VsV%V 2 ± i> (3) 

Here V g = c s /(RQi) is the gravitational drift arising through the magnetic curvature terms, 
V n = c s / (L n fli) is the diamagnetic drift speed, c s is the ion acoustic speed, R is the major 
radius of curvature, Qi is the ion cyclotron frequency, L n is the equilibrium density scale- 
length and /1, and D are the dynamical viscosity and the diffusion coefficient respectively. 
Here Va is the Alfven velocity normalized to the sound velocity (V^f = v^/cf). Thus the 
plasma (3 can be expressed as (3 = The coefficient of resistivity defined as r] s = 

u(c 2 /ujp)(u ci /vl) = (u/uj p )(uj ci /ujp)(c 2 /vl) = v/u) ce is a dimensionaless parameter. Thus, 
r] s = 1.6 x 10 -13 nln(A) / BT^ 2 , here n is the plasma density, B is the magnetic field in c.g.s 



system of units and T e is the electron! temperature in eV, ln(A) is t 



re Coulomb logarithm. 



The model set of Eqs.(^|Hl) has been derived earlier by Kaw (see j^]). In the present work 
we discuss in detail the linear and nonlinear features exhibited by these set of equations. 

Comparing Eqs.Q2J) to our earlier 2d model equations, we see that the extended model 
has additional linear and nonlinear coupling to the magnetic perturbation. The coupling 
coefficient is proportional to V\ (i.e. inversely proportional to (3) as well as to the spatial 
variation in the parallel direction (i.e. to ku). Equation (JHJ) describes the time evolution of 
the magnetic fluctuation and is coupled both to the density and potential fluctuations. In 
terms of basic physics the generalized model has an additional collective degree of freedom, 
namely the shear Alfven modes and finite k» effects bring about a linear and nonlinear 
coupling between them and the RT and drift modes. We also have an additional dissipation 
parameter in the system, namely the resistivity coefficient rj s appearing in the Ohm's law. 
We will discuss the linear properties of the model in greater detail in the next section. 

Note that the two dimensional limit can be obtained by putting d/dz —* for which Eq. (JHJ) 
gets totally decoupled from the equations for n and 0. In this limit, it is easy to show from 
Eq. (JHJ) that the magnetic energy ( ~ / | Vip 2 \ d 3 r ) simply decays away at a rate proportional 
to rjgVj. So in the two dimensional limit the magnetic energy has no role to play in the 
evolution of n and 0. Even when 3d effects are important the electromagnetic effects can be 
negligible. This will happen when the evolution of ip becomes unimportant but the parallel 
current continues to remain finite and provides coupling to finite k z modes in the evolution 
equations for density and potential. Such a limit is possible when 7] s V^V 2 ip » dip/dt (or 
rjsV^k^ >> 7, the growth rate ). The requisite limiting procedure thus consists of letting 
ip — > but letting the parallel current contribution (proportional to on the RHS of 

cq.(JHJ) remain finite. Thus from (JHJ) we have, 

Substituting for V\ip in (JTJ) and (J2J) we obtain, 

d dn - 1 f d 2 1 

ftVV + V g ^ + zxVip- VVV - - l-g^in - <p)\ = M VV (6) 

We will refer to the above two coupled set of evolution equations fEqs. (J5l6|0 as the 3d 
electrostatic model. The set of Eqs. ()5l6|) are considerably simplified in comparison with 
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Eqs. (jl|2|3|) . They describe the coupling between only two variables, viz. density and the 
scalar potential with no electromagnetic effects. However, the influence of the third dimen- 
sion is still present through d/dz for finite values of the parallel component of the wave 
vector fen. Physically, this simplified limit can be understood as follows. The perpendicular 
variation of the Rayleigh Taylor mode produces the charging of magnetic field lines via the 
polarization drift effect. Finite spatial variation in the parallel direction implies that the 
magnetic field line is charged differently at different points thereby promoting the flow of a 
parallel current. The magnetic field associated with this current is responsible for the elec- 
tromagnetic perturbations. However, if the resistivity of the plasma is high (rjsVjkj^ » 7 
the linear growth rate), such a parallel current gets heavily damped, consequently the mag- 
netic field as well as the creation of the rotational electric field dip/dt is negligible (and 
hence the limit i/j — > 0). The perturbations are therefore essentially electrostatic in nature 
in this limit and hence can be considered as an appropriate three dimensional extension 
of the earlier 2d electrostatic model. It provides a simple means of carrying out a direct 
comparison with the 2d results in the presence of finite k\\ and finite resitivity effects. 

In the absence of V g , the gravitational drift, the 3D electrostatic model Eqs. (J5pfi)) reduces 
to the well known Hasegawa Wakatani model [l9| , studied in great detail for the understand- 
ing of electrostatic low frequency plasma turbulence phenomena in three dimensions. Its 
2D variants, obtained by replacing z derivative by a single scalar number has also attracted 
considerable attention |l9l |. 

There is an interesting scaling property displayed by both 3d electrostatic as well as the 
generalized 3d electromagnetic equations which we now wish to highlight. The equations 
remain invariant under the following scaling transformations: 

z^z/a; r] s = a 2 rf s ; ip = V\ = V\ja 2 (7) 

Here a is a scalar scaling factor. These scalings help in establishing equivalence amidst 
a wide class of phenomena for which the parameters j3, r] s and the typical length scales 
along the equilibrium magnetic field direction are related according to the above mentioned 
scaling relations. Note that the transformation leaves k z VA, ^Va and rjsV^ invariant. It 
also leaves the total energy (as well as each of the individual components of energy, namely 
pressure, kinetic and the magnetic energy) as invariant. Although the field ip gets scaled, 
yet the magnetic energy (normalized to plasma thermal energy viz. b 2 /8irnT ) which is 
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f(Vip) 2 d 3 r/(2f3 J d 3 r) = V\ f(Vip) 2 d 3 r/(2 J d 3 r) in our normalizations, remains invariant. 
The scaling relationship helps in carrying out simulation for a convenient choice of the set 
of parameters j3, t] s and L z (the box length along the z direction which defines the typical 
size of the excitation scales along z), which can later be related to the realistic set of values 
using the scaling coefficient a. 

We will explore the nonlinear states of the 3D electrostatic model set Eqs. (l5l6J) as well 
as the full generalized electromagnetic set given by Eqs. (jll2l3|) in sections IV and V after 
discussing their linear properties in the next section. 



III. LINEAR ANALYSIS 



The coupled set of equations (I1I2I3I) can be linearized and fourier analyzed to obtain the 
following dispersion relation 

- ik\uo 3 + {{D + fi + rj s Vl)ki + iklky(V g + V n )} uo 2 
+ {i{Dq a VX + Dfi + 7i.V%n)kl + %k 2 z V 2 A k\{\ + k\)}uo 

+ {-kiky [{r) s Vl + fi)V g + (D + fi)V n ] + tk 2 y (V 2 - VgV n - k 2 ± V g V n )}u 
- DrjsV^kl — (D + iik\)k 2 z V 2 A k\ - ifik b ± k y ( Vs V^V g + DV n ) + ik 2 ± k 2 z V%k y (V g - V n ) 

- r) s Vlklk 2 y V g {V g - V n ) + fik\k 2 y V n V g - ik\V 2 V n + ik 3 y V g V 2 = (8) 

The above dispersion relation contains three basic modes, namely, the drift wave, the 



Rayleigh- Taylor mode and the shear-Alfven wave. This can be seen quite easily by set- 
ting all the dissipative coefficients to be zero (i.e. D = fi = rj s = 0) and by rearranging 
Eq.(JHJ) in the following form. 



kyV g uj + iv g (V n -V a )} = k 2 z V 2 {l + k] 



ky{V g - V n ) 



{u - k y V n ) { u, 2 - kyVgU + £ V g ( V n - V„) \ = k 2 z Vi(l + ki){u+ ^ 9 - k2 J I (9) 

For V g — V n — (i.e. in the absence of magnetic curvature and density gradients) Eq.® 
gives the kinetic Alfven wave dispersion relation. 



OJ 



k 2 V 2 (l + k 2 ± ) (10) 



When V n is finite, V g = and V\ — > oo (a low (5 plasma) we have upon dividing Eq.© by 



V 2 



A- 



k V 

V n I -i -i \ 
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which is the drift wave dispersion relation. For this case the electrons have a Boltzmann 
distribution, i.e. the wave time scales are in the regime of u / k z Vth, e < 1- The two dimensional 
electrostatic Rayleigh Taylor growth rate can be recovered by putting k 2 2 V\ — > in J0J). 



For this mode the electrons act like a two dimensional fluid, under the condition of 
u/k z v t h. e > 1- This is the only mode which is unstable (has a finite growth rate) in the 
nondissipative limit. At a finite value of k z VA this unstable mode gets coupled to the stable 
Alfven branch. We show in Fig.l the variation of the real and imaginary parts of u as a 
function of k z VA- The value of k x = k y = 0.1 has been chosen for the plot. The other 
parameters are V g = 0.036, V n = 0.8, D = /i = r] s = 0. For small values of k z VA we observe 
that the three roots of the cubic equation are essentially obtained by putting the right hand 
side of Eq.(jnj) to be zero. This gives rise to one real root uo = k y V n (arising due to the 
balance between the parallel electric field and the equilibrium pressure variation along the 
bent magnetic field lines ), and two complex roots of the Rayleigh Taylor mode. For the 
parameters k « 1 (scale lengths longer than p s ) we have (cu = u> r + iujj) 



Clearly, for a choice of V g « V n the real part of the frequency is much smaller than the 
growth rate i.e. uj r « Ui and is negligible as the plot of Fig.l in the regime of small k z VA 
shows. As k z VA increases the right hand side of Eq.fjHJ) cannot be ignored. Figure 1 shows 
that with increasing k z VA, \ &i | decreases and goes to zero at k z VA = k zc VA- The point 
k z = k zc is in fact the point of exchange of instability, as u r w remains close to zero upto 
this point for the Rayleigh Taylor branch and becomes finite for values of k z VA > k zc VA- 
The expression for k zc can thus be determined from Eq.Q by substituting uj = 0. This gives 



The critical wavenumber k ZCl beyond which the growth rate vanishes, thus increases with 
increasing values of V n , V g and k y but decreases with increasing values of Va and k x . For 
the values chosen for these parameters in Fig.l we have k zc VA = 0.12. The limit on k z 




(12) 






(13) 
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for instability is in fact identical to the threshold condition on plasma (3 encountered in 
the context of ideal Magnetohydrodynamic (MHD) ballooning modes in toroidal devices 
like tokamaks. In the ballooning mode case, one generally seeks the critical value of of 
plasma beta ( (3 = /3 C ), beyond which the instability sets in for a fixed value of the parallel 
wavenumber. Here, on the other hand, we have fixed the value of f3 and are seeking the 
threshold condition on k z the parallel wavenumber below which the instability exists. This 
has been done keeping in view the identification of the linearly unstable modes for the three 
dimensional simulations, where a range of k z are present and (3 = 1/V| would be taken 
as a parameter. However, the expression for (3 C encountered in the context of ballooning 
modes can be recovered from Eq.(JT3J) by substituting k zc = 1/qR, k y ~ k±, 1/V| = /3 C , 
V n V g = l/RL n = 1/Ra, (as the density gradient scale length L n can be taken typically to be 
of the order of minor radius a). These substitutions in Eq. (fT^|) then lead to the well known 
expression for the critical value of plasma beta as (3 C = a/q 2 R = e/q 2 , where e is the aspect 
ratio. 

In the region k z < k zc the curve u>i vs. k z in Fig.l typically seems to have an elliptical 
shape. This can be understood as follows, we know that in this region to r — > 0, hence 
uj 2 = —ujf. Considering the regions where Ui < k y V n and using V g « V n we write the 
dispersion relation as 

^ + k 2 z V 2 -^V g V n = 0- i.e. u 2 + {k 2 z - k 2 zc )V 2 = 0. (14) 

which is an equation of an ellipse. Interestingly, the dispersion relation simplifies to similar 
elliptic form as of Eq. fpHj) for uji > k y V n , with the only modification that in this case the 
term k 2 z V\ is multiplied by the factor of (1 + k\). 

There are certain other features exhibited by the plot of Fig.l. The only root which has 
the finite real part at k z = is u r = k y V n , and arises from the decoupled ip equation in this 
limit; at higher k z it approaches the stable shear Alfven branch with uj r = k z VA- This mode 
remains stable throughout, the imaginary part of this particular mode remains = in the 
entire k z domain. The real part of the both RT branches are zero for k z < k zc ; however as 
k z is increased beyond k zc , one amongst them asymptotes towards the drift wave dispersion 
relation u r = k y V n /(l + k 2 ^) and the other approaches the complementary branch of the 
shear Alfven mode i.e. u = —k z VA- 

We now investigate the effect of dissipation on the frequency as well as on the growth rate 
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of the three modes. In Fig. 2 we have plotted the growth rate and the real frequency with 
1z z Va when the dissipative coefficients /i (subplot (a) and (b)) and r] s (subplot (c) and (d)) are 
separately taken to be finite. A comparison with the corresponding non dissipative case (all 
other parameters being identical) of Fig.l, clearly shows that a finite value of [i does not alter 
the real frequency and the threshold condition on k z VA significantly. It, however, causes an 
overall reduction in growth rate. An entirely different and interesting trend is observed with 
respect to the dissipative coefficient r] s . As the value of r] s is increased the real frequency 
instead of suddenly acquiring a finite value beyond k z VA = k zc VA, gradually starts deviating 
from zero even before k z = k zc to finally asymptote towards the drift k y V n /(l + k'jj and 
the shear Alfven branch —k z VA at large enough k z VA- The growth rate too unlike the case 
of rj s = continues to remain finite even for k z > k zc i.e. the unstable domain of the k z 
space gets widened. This happens essentially because of the presence of resistivity driven 
modes; viz. the resistive 'g' and the resistive drift mode. The expenditure of energy in 
causing the field line bending along the parallel direction (for u/k z v t h e > 1 modes ) leads to 
the stabilization of the RT instability for finite k z . The field lines bend due to the parallel 
component of the current. In the presence of resistivity the parallel currents gets damped. 
This leads to the recovery of instability in the resistive time scales, and is the physical basis 
of the excitation of the resistive interchange mode. The resistive drift wave arises in the 
regime u < k z v t he where the n, (p relation wants to be Boltzmann like but acquired deviations 
because of finite resistivity effects; the phase difference between n and vph permits energy 
exchange between the waves and the fluid and leads to the resistive drift wave instability. 

In Fig. 3 we depict the features of resistivity driven branch in more detail. It shows the 
plot of real (solid lines) and imaginary (dashed lines) part of the frequency for this branch. 
The ideal results t] s = are plotted in subplot(a) (V g = 0.0364) and subplot(c) (V g = 0.01) to 
be compared with the finite resistivity r] s V2 = 1.6 plots of subplots (b) and (d) respectively. 
As mentioned earlier a finite growth rate beyond k z = k zc is because of the two resistivity 
driven modes, namely the resistive interchange mode and the resistive drift wave. The region 
where Ui » uo r has the characteristics of the resistive interchange mode; in the opposite 
limit u r » uj,i it is the resistive drift wave j^J. The figure also shows that as we reduce 
the value of V g the resistive drift regime gets broadened, (u r asymptotes to the drift wave 
frequency at a lower k z ). 

In Fig. 4 we have isolated the growth rate of the resistive drift mode by choosing V g = 
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0. One can see from the plots that the growth rate due to this mode vanishes for k z = 
0. An analytical expression for the growth rate of this particular mode can be obtained 
perturbatively in the small rj s limit as we show below. We look at the dispersion relation 
of Eq.fjSJ) in the limit of finite r) s . We put D = fi = for simplification. The dispersion 
relation in this case is similar to Eq.® but with an additional rj s dependent term, which is 
responsible for instability. 

(u - k y v n + i Vs klvi) jc 2 - k y v g u + &v 9 (v n -v g )} = e z vj(i + kl) L + ky ^lJj ] } 

X (15) 

The effect of 7] s on the drift wave can be seen by taking the limit V\ — >• oo and V g = in 
Eq.(HSD. 

k y V n . Tj s k^ 2 . . 

(16) 

Considering the term on right hand side as a small correction (possible when rj s is small and 
k z is large) the dependence of rj s on uj can be obtained iteratively as 

U= (l + fc 2) +% k 2 (l + fc 2)3 

which shows the resistive destabilization of the drift wave branch. The expression obtained 
above shows that the growth rate reduces as k z is increased. The small k z limit can be 
captured by ignoring uj in the quadratic dispersion relation of Eq. (fTH|) giving 

—mm 

This shows that the growth rate increases with k z indicating 7 vs. k z must pass through a 
maximum, as seen in Fig. 4. 

We next look at the linear properties of the simplified three dimensional electrostatic 
model equations (Eq. (|5l6j )). This model basically contains the Rayleigh Taylor mode and 
the drift wave mode, but the coupling to the shear- Alfven branch is absent. The dispersion 
relation for this set is 

- 2 + »{-kyV g + *f + - V 9 ) - (K -V a ) = (19) 

Clearly, one recovers the two dimensional RT growth rate expression in the k\ = or 
r] s — * 00 limit. The first order perturbative corrections show damping due to r) s and add a 
perturbative correction to the real frequency of this mode. The drift wave dispersion relation 
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can also be recovered provided k z is finite. Taking V g = and considering rj s to be small 
the dominant balance is between the two k z dependent terms in f!19|) from which we obtain 
the standard drift wave dispersion relation of u = k y V n /(l + k 2 ). Retaining the first order 
correction due to the remaining rj s dependent term we obtain the same expression for the 
resistive destabilization as that of Eq.(|17|). 

The role of r] s on the drift wave branch for the general non perturbative case can be gleaned 
from the plot of Fig. 4 for both the electromagnetic and the simplified electrostatic models. 
The upper subplot for r] s V^ = 1.6 clearly shows that the two growth rates obtained from the 
two models differ at small k z VA but agree well at large values of k z VA- At small k z VA the 
growth rate obtained from the fully electromagnetic dispersion relation is found to be smaller 
than the electrostatic case, which basically shows that the addition of electromagnetic effects 
cause stabilization. The two growth rates, however, show a similar trend in both cases, 
namely, they first increase with k z and reach a maximum value and then fall off with k z . 
The decreasing trend with k z is captured by the perturbative expression of Eq. fTTj) . which 
is valid when ri s /k z is small. Furthermore, the lower subplot of the figure also shows that 
for large values of T]^/\ (= 16) the agreement between the two models is excellent over the 
entire range of k z VA- This is a further evidence of the fact that the approximations used 
in the derivation of the simplified electrostatic model are very accurate in the limit of large 
resistivity. Though the other two modes are essentially damped by the finite value of r) s , 
the resistive destabilization of the drift mode can enhance the growth rate of the coupled 
system. 

To summarize, linear analysis of the system shows that inclusion of the third dimension 
introduces additional unstable modes and the presence of electromagnetic effects brings 
about a coupling to shear Alfven modes. It is then of interest to understand the nonlinear 
evolution characteristics of these modes and their evolution into possible saturated nonlinear 
states. The simplified electrostatic model can be useful in isolating the physics due to finite 
k z (parallel variation ) from the electromagnetic characteristics. In the next two sections 
we present the nonlinear evolution studies of these two models with the help of numerical 
simulations. 
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IV. NONLINEAR SIMULATION RESULTS FOR THE 3D-ELECTROSTATIC 
MODEL 



We begin by presenting and discussing the results of numerical simulation studies of the 



). The equations are evolved with 
Jl. Most of the studies have been 



simplified 3d-electrostatic model represented by Eqs. 
the help of a fully dealiased pseudospectral scheme 
carried out with a resolution of 64 x 64 x 64 fourier modes in the three directions. Some 
test studies were also carried out with a lower resolution of 64 x 64 x 16 modes. 

We first investigate the question of the existence of saturated nonlinear states. In our 
earlier two dimensional studies it was shown that, even in the absence of any boundary 
or initial condition related anisotropy, there are two distinct symmetry states to which the 
system goes in the nonlinear state. The system ultimately forms growing streamer structures 
which are radially (along x in the slab description) elongated or to poloidally (y direction) 
symmetric saturated zonal patterns. For a fixed value of the driving parameters V n and 
V g , the condensation to these symmetry states was governed by the value of dissipative 
coefficients D and \x. We now investigate the role of three dimensional perturbations on the 
development of these symmetry patterns. We choose the perpendicular box dimensions as 
L x = L y = 20np s ; so that there is no boundary related anisotropy in the perpendicular plane. 
L z is chosen as 1.25 x 10 5 p s ; this is to concentrate on low frequency modes of interest which 
are extended along the field lines and have k z /k± << 1. We choose the parametric regimes 
close to the 2d case in order to study the effects of parallel scales k z and the resistivity 
parameter r] s on the formation of nonlinear states. We chose D = \x = 0.1, V n = 0.8, 
V g = 0.036, a parametric region that corresponds to saturated zonal patterns for the two 
dimensional case and have made several simulation runs for different values of r] s ( < 10~ 5 ) 
and starting with small amplitude random initial perturbations in n and (p. The chosen 
values of V n and V g also correspond to those of the currentless toroidal device BETA on 
which several experiments on the MCD-RT have been carried out. For this machine the 
typical value of rj s = v/uj ce ~ 10 -5 and the parallel scale lengths are typically of the order 
of L z « 125Mts. ( L z = 2imR with the major radii R m 45cm. and n w 30 — 50 is the 
observed toroidal winding number of the magnetic field lines in this machine). Our results 
are applicable to other values of r] s and k z through the scaling arguments described with 
Eq.(0) and may also be applied to the spread F - region of the ionosphere, ELM region of 
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tokamaks etc. (for more details please see the last section). 

A typical 3d representation of the initial random potential <p structure is shown in the 
form of a slice plot in Fig. 5. These slice plots basically show the pattern of a particular 
field variable with the help of color (in color plots) or through shading (in gray plots) on 
various two dimensional slice planes of the three dimensional space. The appropriate slicing 
helps in the three dimensional visualization of the field pattern. For instance in the plot of 
Fig. 5 (and also in all the subsequent slice plots presented in this paper) we have chosen five 
different slices of the three dimensional volume of L x x L y x L z . The different slices are the 
x vs. y plane at z = and z = L z /2; y vs. z plane at x — L x /2 and L x ; and x vs. z plane 
at y = 0. 

The initial density n field has a similiar structure containing many short scales in a 
random configuration. For the small amplitudes chosen initially the evolution is primarily 
governed by the linear terms. Thus, during the initial linear phase modes having maximal 
growth rate acquire the largest amplitudes and dominate the spectrum. Since the maximally 
growing modes are the RT modes that have large k y but small k x , we can expect to see, in the 
early linear stages of development, the appearance of structures which are elongated along 
the x direction. We see clear evidence of such structures in the slice plot for the potential 
at t = 30 in Fig.5. As the amplitudes grow, the modes start interacting and one expects the 
power to get transferred to linearly stable modes as well. In the present case we see such 
a phenenomenon too and find the power in potential tp field nonlinearly cascading towards 
long scales. Such a cascade towards long scales is an intrinsic property of the polarization 
drift nonlinearity that is present in the evolution equation of tp. 

We have carried out simulations both for large and small resistivity parameters t] s , and 
observe distinct difference in the two regimes. For large values of r) s , viz. 10~ 5 , nonlinear 
saturated zonal symmetry patterns in potential </? field are seen to form (see Fig.5). However, 
when T] s is small = 1.26 x 10~ 6 there is no saturation and growing streamer patterns are 
observed (Fig.5). Thus in three dimensions, as the available phase space of the modes get 
enhanced with the addition of finite k z modes, the resistivity parameter 7] s along with D 
and /i determine the symmetry pattern of the potential structure in the nonlinear stage. 
Fig. 6 shows the simulation cases in the parametric space of D vs. l/r] s for which saturated 
states were achieved (by circles) and those for which only growing streamer patterns were 
observed by + (plus) signs. This trend is consistent with the fact that at large values of rj s 
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the contribution of the additional linear and nonlinear terms become small in comparison 
with other terms and the set of equations tend to reduce to the previous two dimensional 
equations. This has been quantitatively illustrated in Fig. 7, which shows the plot of the 
ratio of growth rate 7 with k^/r] s as a function of k z for the two values of r) s yielding 
saturated zonal (dots) and growing streamers (+ sign). The range of k z shows the permissible 
parallel wavenumbers of the simulation after aliasing. It is clear from the plot that when 
rj s = 1.26 x 10 -6 , there are parallel scales for which the ratio drops below unity (signifying the 
dominance of the extra three dimensional terms in the evolution equation and consequently 
the dynamics being altered significantly. On the other hand for large values of t] s = 10 -5 
there are no parallel scales for which the additional terms dominate, the dynamics thus is 
close to the two dimensional scenario yielding saturated structures. 

However, there are a few interesting differences in the composition of the final saturated 
state for the two and three dimensional cases even though the addtional terms are merely 
small perturbative corrections for such numerical runs. In the two dimensional simulations 
the density field was observed to get slaved to tp and it too displayed the formation of long 
scale structures with two distinct symmetries. For the three dimensional runs however, we 
see, from the plot of Fig. 8, that the density field continues to be dominated by power in 
the short scales. The scatter plots between the density and the potential fields (see Fig. 9) 
also does not show any evidence of functional relationships developing between density and 
potential fields. The vorticity V 2 y? too, unlike the previous case, has significant power in 
short scales (Fig. 8) and does not form any functional relationship with the potential tp field 
(see Fig. 9). 

The non slaving of the density field can be understood by realising that due to parallel 
variations, additional modes (drift waves etc.) having very different linear mode relationship 
amidst the two fields compared to the RT modes get excited (e.g. typically rik ~ Pk for drift 
waves, whereas for RT mode the density and potential fields are essentially out of phase ). 
This may hinder the slaving process of density to the potential. Moreover, two dimensional 
set of equations conserve the following non dissipative integral invariant 



which clearly shows an establishment of integral relationship between the density and the 
potential fields. Such a integral relationship is consistent with the possibility of the n field 




(20) 
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getting slaved to <p. The incorporation of 3d effects, however, rules out any integral constraint 
on the two fields. We have in this case instead: 

(21) 

The dissipation parameter 7] s has been retained, as the 3d effects here arise solely from r) s 
dependent term. Clearly, since the two fields do not satisfy any integral relationship in the 
3d case, it is not easy for the density field to get slaved to tp. It thus seems that the non 
existence of any integral constraint and the absence of any functional relationship between 
n and ip along with the fact that the nonlinear evolution of the density field is governed 
directly by the convective nonlinearity viz. z x S7<p ■ Vn which cascades power towards short 
scales, leads to the predominance of power in short scale fluctuations in n. The polarization 
nonlinearity influences the n evolution only indirectly through <p. 

It should be noted that in the above runs with finite r) s , the linear phase has two unstable 
modes - the RT mode and the drift mode which is made unstable by the resistivity. In 
order to understand the role of the unstable drift wave it is possible to isolate its behaviour 
by artificially turning off the RT mode. We have carried out such an investigation by 
setting V g = and looking at the nonlinear evolution of the drift modes. The value of the 
resistivity parameter was chosen to be rj s = 10~ 5 . Figure 10 shows a comparison of the 
growth and evolution of the total energy, the zonal and streamer powers for RT mode and 
the resistive drift wave. The slower linear rise can be understood from the lower growth rate 
of the resistivity driven drift wave in comparison to the growth rate of the RT mode. It is 
interesting to observe that in the final saturated regime of the resistive drift case, there is 
no dominance of power in the zonal mode as observed in the context of RT. This leads to a 
characteristic mixed flow pattern in which one cannot clearly distinguish between the zonal 
and streamer symmetries. Such a saturated state of the potential fluctuation at t = 450 is 
shown as a slice plot in Fig.ll. It is also interesting to observe that for the resistive drift 
wave the energy level in finite k z modes is an order of magnitude higher than that in the 
k z = modes. 

In all of the above simulations we have avoided introducing any perpendicular anisotropy 
associated with boundary and initial conditions. We restricted ourselves to those simulations 
for which the perpendicular aspect ratio of the simulation box was unity i.e. L x = L y . The 
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simulation volume L x x L y x L z basically represents a small region of the entire plasma. By 
simulating over a small region one hopes to identify and understand the basic features of 
turbulent excitations. The simulation, however, gets constrained by the choice of box sizes 
in a few ways; the longest scale length along a particular direction is determined by the box 
length along that direction. In some cases as we would show below, a natural process of 
power cascade towards a long scale asymmetric mode can get inhibited by a certain choice 
of the aspect ratio of the simulation box size. The unlimited growth of streamers observed 
by us earlier for the choice of aspect ratio unity is one such example. It is well known that 
a shear flow excites Kelvin - Helmoltz (KH) instability; however, in our simulations carried 
out with with L x = L y we observe no development of secondary KH instability which could 
prevent the unlimited growth of streamers in the parameter domain indicated by the + sign 
of the plot in Fig.6. This happens because the essential condition for the excitation of KH 
instability on streamer shear flow (with shear scale length of L y ) can never be met within the 
restriction of square box size. The KH instability can be excited only when the perturbation 
scales (in the orthogonal direction) are longer than the background shear scale length. With 
L x = L y there can be no such mode to support such a secondary destabilization process of 
streamers. By relaxing the constraint of L x — L y , and choosing instead L x > L y (we chose 
L x = 4:L y ) we observe that the unlimited and unphysical growth of streamers is prevented. 
In Fig. 12 we show a comparison of the evolution of total energy / / (n 2 + (Vip) 2 )dxdy for the 
two cases, viz. L x = L y (solid lines) and L x = 4L y (dotted lines). The other parameters for 
the two cases are identical. The plot clearly shows that the when the box dimensions are 
identical the energy grows indefinitely and there is no saturation. On the other hand when 
L x is chosen to be longer than L y the energy saturates. This has important implications on 
transport. It shows clearly that simulations with those parameters which require L x > L y 
for saturation; excites flows with radial scale length longer than the ones for which the 
instability saturates for L x = L y . Since the radial decorrelation step size is essentially 
determined by the radial scale length of the structures, it implies that the transport will be 
high for the parametric regime (+ sign of Fig.6) which do not saturate for a square box size 
(aspect ratio unity). 

To summarize, in this section we have shown that similar to the 2d simulations the 3d 
electrostatic model too is capable of supporting nonlinear saturated states that are domi- 
nated by long scale zonal flow patterns. Here too, even in the absence of any boundary or 
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initial condition related bias, the nonlinear evolution lead to the condensation towards either 
saturated zonal flow patterns or growing streamer formations depending on the strength of 
the various dissipation parameters. The model has an additional dissipation parameter in 
the resistivity coefficient which plays a special and distinctly different role than the other 
two dissipation parameters in the selection of the final stage. The density field in the 3d sat- 
urated cases show distinct features of short scale dominance and non slaving to the </? field; 
which is distinctly different from the 2d results. Furthermore, the resistivity parameter also 
leads to the excitation of additional instabilities e.g. resistive interchange and the resistive 
drift waves. Simulation studies on the unstable resistive drift wave were also carried out 
which show new variety of nonlinear state in which zonal and streamer powers were com- 
parable. It was also shown that the parameter regimes for which one obtains unsaturated 
streamer patterns could be stabilized by increasing the aspect ratio L x /L y from unity. This 
basically enables the secondary KH destabilization of the streamer patterns. It was shown 
that such cases would be responsible for higher transport. 

V. SIMULATION OF THE 3D ELECTROMAGNETIC MODEL 

We now turn to the generalized three dimensional model represented by Eqs.(JI]|HJ) and 
discuss its numerical solutions. In comparison to the 3d electrostatic case we now have 
an additional field variable ip (the magnetic vector potential), whose temporal evolution is 
governed by Eq. and which provides additional coupling (linear and nonlinear both) terms 
in the evolution equations of the density and potential variables. The coupling coefficient is 
proportional to V\ (i.e. inversely proportional to (5) as well as to the spatial variation in the 
parallel direction (i.e. to fen). As mentioned before the influence of electromagnetic effects 
will be felt when dip/dt becomes comparable to rjgVj'V 2 '^. The order of typical perpendicular 
wavenumbers ranging from 0.1 to unity for our simulations, it implies a direct comparison 
of r] s Vjl with the growth rate. When rj s V^ » 7, ip is essentially damped, electromagnetic 
efects are weak and energy in the magnetic field is typically small. This is the limit where 
the generalized 3d electromagnetic model is expected to reduce to the 3d electrostatic model 
discussed in detail in the last section. The parameter regime of rj s V\ « 7 to t] s Vj of the 
order of 7 is thus of interest for studying the influence of electomagnetic effects on the 
dynamics. 
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In this section we present results for the case of rj s Vj[ = 1. For this value of r] s V2 
we expect the electromagnetic effects to become significant enough so as to influence the 
dynamics. Furthermore, this choice is also motivated by the fact that it is the parameter 
regime relevant for the currentless toroidal BETA machine at IPR. A detailed parametric 
simulation study for t] s V2 ranging from << 7 to >> 7 for the 3d electromagnetic set of 
equations are, however, underway and will be presented elsewhere. 

We look for saturated states in a square box geometry with (L x = L y ). The other 
parameters are V g = 0.036, V n = 0.8, D = fi = 0.1 and 2ttVa/L z = 10~ 3 . The plot in 
Fig. 13 shows the evolution of the total and the magnetic energy (solid and dashed lines 
respectively). The energy is seen to saturate after an initial exponential growth. The 
saturated magnetic energy is an order of magnitude smaller than the total energy for this 
set of parameter values. Thus the electromagnetic fluctuations are not dominant in this 
particular smulation. However, we present a comparison with the electrostatic case to show 
that even in this case, though the amplitude of electromagnetic fluctuations are low, their 
effect on the nonlinear saturated state is significant. The linear growth of energy for both 
electromagnetic and electrostatic cases are similar. This is because the value of the maximum 
growth rate in both the cases are identical. However, in the nonlinear regime the electrostatic 
total energy is smaller due to a mild decay during the later phase. The lower subplot of the 
same figure shows the evolution of intensity of zonal and the streamer modes. The zonal 
intensity for both electrostatic and the electromagnetic cases are at identical level. The 
intensity of the streamer mode, however, in the electrostatic case is considerably smaller 
and exhibits a mild decay similar to what is observed for the total energy. It is interesting 
to note that in the electrostatic case even though the total energy in turbulent fluctuations 
are small compared to the electromagnetic case, the zonal intensity is identical, i.e. as 
a relative fraction, zonal flow intensity is stronger in the electrostatic case. This clearly 
implies that it is much easier to generate zonal flows in the electrostatic case, confirming 
the prevalent lore that that electromagnetic effects inhibit the zonal flow generation. The 
streamer intensity in electromagnetic case being higher also confirms that the stabilization 
of the instability in the presence of electromagnetic effects becomes difficult. 

In Fig. 14 we depict the slice plots for the three fields viz. n, ip and ip in both the linear as 
well as the nonlinear regimes. In both linear and nonlinear regimes we observe considerable 
structure in all the three fields in the z direction, indicating that energy in the finite k z modes 
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are significant. Unlike the clear symmetry breaking nonlinear stages (growing streamers or 
saturated zonal depending on the parameter regimes in the D, \x and rj s space) of the 3d 
electrostatic RT modes, in the electromagnetic case the flow structures cannot be distinctly 
classified in zonal and streamer patterns. This is also evident from the evolution of zonal 
and streamer intensities in the plot of Fig. 13. Furthermore, the n and ip fields show a 
predominance of short structures in comparison to the potential (f field. 

Our simulations also show that similar to 3d electrostatic case, here too the density re- 
mains an independent field throughout the evolution and does not get slaved to the potential 
field. In the 3d electrostatic context, this was attributed to the the presence of a variety of 
additional modes arising by permitting the three dimensional variation in the system and 
also to the loss of integral invariant for the 3d equations. In the electromagnetic case where 
there is an increase in the variety of linear modes, (shown in detail in the third section ) it 
is even more difficult for the two fields to develop any functional relationships. 

The electromagnetic studies of this section reveal that the features, which were earlier 
(with the help of electrostatic studies), attributed to the three dimensionality of the sys- 
tem are present in these simulations also, the sytem being three dimensional here as well. 
However, additionally we observe that even the presence of a weak electromagnetic energy 
considerably opposes the process of nonlinear stabilization. 



VI. SUMMARY AND CONCLUSIONS 



In this work we have studied the extension of a previous two dimensional nonlinear model 
evolution equations 0,1^1 for the magnetic curvature driven Rayleigh Taylor instability to 
three dimensional perturbations. The extended model also incorporates coupling to elec- 
tromagnetic fluctuations associated with the magnetic field line bending terms through the 
parallel component of the Ohm's law. The objective of the present work has been to inves- 
tigate the influence of three dimensionality and the electromagnetic effects on the nonlinear 
state. It was shown that the effects due to three dimensionality can be isolated by consid- 
ering a simplified 3d electrostatic limit. Such a limit is valid when the typical growth rates 
are much smaller than the parameter r\ s k\V\. Here rj s is the resistivity parameter defined 
earlier in the text, f3 = 1/V| represents the plasma beta and k±, the typical perpendicular 
scales. 
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Studies on 3d electrostatic model show that the D — p phase space, which in earlier 2d 
studies jl| governed the symmetry of nonlinear flow patterns (viz. the transport inhibiting 
saturated zonal flows versus the transport enhancing growing streamers) gets extended by the 
inclusion of a third dissipative parameter viz. r) s . A comprehensive parametric study reveals 
that the resistivity parameter, has an entirely different role in pattern selection process. 
While smaller values of both D and p form zonals and their larger values streamers, it 
is the opposite for rj s . In 2d simulations it was observed that the density field ultimately 
develops a functional relationship with potential. The power cascade towards long scale for 
the potential in 2d case ultimately also forces the density to acquire long scale structures. It 
was observed in our current 3d simulations that the density was not in any way constrained 
to follow the potential field. Hence, in the 3d case the density field continues to have power in 
short scale fluctuations. Another feature of the 3d model, in contrast to 2d, is the existence 
of additional resistivity driven modes. Nonlinear studies on resistively destabilized drift 
wave (with RT growth switched off) yielded a novel variety of saturated states which have 
neither zonal nor streamer symmetries; instead the intensities of both zonal and streamer 
flows were at comparable level. 

The electromagnetic effects on RT were studied by simulating the fully 3d electromag- 
netic set of equations. At the moment we have carried out investigations only for those 
parameters for which the electromagnetic effects are weak and the magnetic energy is an 
order of magnitude smaller than the total energy. This was achieved by choosing r] s Vj = 1, 
which is still larger than the maximum growth rate 7 maa: = 0.16. These simulations clearly 
show that even though the electromagnetic energy is weak, its presence inhibits the zonal 
formation leading to a relatively higher amplitude of streamers. Thus, as expected, the 
presence of electromagnetic effects hinders the saturation process. 

Curvature or gravity driven RT modes are important in many magnetized plasma prob- 
lems. They were extensively studied experimentally in the toroidal currentless plasma ma- 
chine BETA (21]]. As stated in the text, the parameters of BETA were such (V n = 0.8, 
V g = 0.036, L z = 1.25 x 10 5 p s , r] s ~ 10~ 5 ) that the results of present study are directly appli- 
cable. In the context of tokamaks, curvature driven ballooning modes (with or without resis- 
tivity effects ) are relevant to the core region, as well as the edge region (especially when the 
plasma ehibits the edge localized modes called ELMS). The typical parameters in the core re- 
gion are (rj s ~ 10~ 7 , V% = l//3= 10, implying that t] 8 V\ ~ 10" 6 , L z /p s = 2nqR/p s = 2xl0 4 , 
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V n = p s /L n « 0.1 and V g = p s jR ~ 10 2 ) and those in the edge region are (r] s ~ 2.5 x 10 5 , 
VI = 1/(3 = 10 3 - 10 4 , so that ri s V\ ~ 1, L 2 /p s = 27rgi?/p s = 2 x 10 4 , K = p s /L n « 0.1 
and = p s /i? ~ 10~ 3 ). Thus the present study may be applicable to the tokamak 
edge problem, although the driving instability in the simulations is somewhat stronger be- 
cause of the higher value chosen for V n . Gravity driven RT modes are also relevant to the 
spread - F region of the ionosphere. The typical parameters in this region are (j] s = 10 -5 , 
Vl = 1/(3 « 10 6 , r) S Vl ~ 10, L zPs = 10 5 , v n = p s /L n « 2.5 x lO" 3 and V g = 10" 4 ). Thus 
the results of our studies are applicable to this problem also. 

It should be pointed out here that several three dimensional numerical studies on the 



electrostatic 
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23j and electromagnetic |2J,|25[ nonlinear equations describing the coupling 



of drift ballooning modes have been carried out in the last decade or so. Such studies have 
employed a realistic three dimensional model for tokamaks with effects due to magnetic shear 
and parallel flows. These studies contain extensive details of simulation studies on transport 
in tokamaks. Our objective in contrast has been to delineate the parameter regime for the 
formation of transport inhibiting and transport enhancing structures and identification of 
the rudimentary physics with the choice of a simplified model ( that of drift - Rayliegh 
Taylor coupling, in the absence of both parallel flow and magnetic shear). 

Finally, we make some remarks on the further exploration of the present work that we 
are currently pursuing. The complete parametric study for the fully 3d electromagnetic case 
has not been presented here. Future studies will explore other regions of parameter space 
(especially lower values of i] s V\ where electromagnetic effects become more important). Thus 
future investigations will be carried out to understand the nonlinear stabilization process by 
the zonal flow formation as the the parameters r\^/\ and f3 = 1/V^f are varied. We also want 
to include the finite Tj effects, bacuse finite ion Larmor radius stabilization is an important 
linear mechanism of the stabilization of the curvature driven instabilities. 
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FIGURE CAPTIONS 



Fig.l Plot of the real frequency (u r ) and the growth rate u>i vs. k z VA for the electromagnetic 
dispersion relation of Eq.(|5j) in the non-dissipative limit i.e. for D = fi = r] s = 0. The 
other parameters are V g = 0.036, V n = 0.8. The perpendicular scales are k x — ky — 0.1. 

Fig. 2 Plot of real frequency and the growth rate from the dispersion relation of Eq.fjHJ) as 
a function of k z VA when fi = 3 is finite (subplot(a) and (b)) and when r\<y\ = 1.6 is 
finite (subplot (c) and (d)). The other parameters are same as that of Fig.l. 

Fig. 3 The plot of real (solid lines) and imaginary (dashed lines) part of the frequency for 
the resistivity driven branch. The ideal results for rj s = are plotted in subplot (a) 
(V g = 0.0364) and subplot(c) (V g = 0.01) for the purpose of comparison with the finite 
resistivity r\ s V\ = 1.6 plots of subplots (b) and (d) respectively. The other parameters 
are the same as that of Fig.l. 

Fig. 4 The two subplots show growth rate vs. k z VA for resistive destabilized drift wave 
V g = 0; (the other parameters being V n = 0.8, , k x — k y — 0.1, D — /i — 0) from the 
fully 3d electromagnetic dispersion relation of Eq.fjHJ) (solid lines, here ) and for the 
simplified 3d electrostatic dispersion relation of Eq.([19j) (circles). The growth rates for 
the two models differ at small k z VA in the upper subplot for which resistivity parameter 
t) s Va = 1-6. They yield almost identical growth rates for large t] s Va = 16 as shown in 
the lower subplot. 

Fig. 5 The slice plots (described in text) for potential <p field to visualize its three dimensional 
structure at various times for the 3d electrostatic evolution equations (0 El)- The 
parameters for this case are V n = 0.8, V g = 0.036, D = fx — 0.1. The simulation 
box sizes are L x = L y = 20np s ,L z = 1.25 x 10 5 p s The first three plots (as indicated 
on the top) are for r] s = 10 -5 at t — 0, 30an<il50 and show the formation of zonal y 
symmetric pattern in the final nonlinear state at t — 150. The fourth plot corresponds 
to a different run for which r] s is = 1.26 x 10~ 6 , and shows formation of radially (along 
x) extended streamer structure. 
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Fig. 6 The parametric regime of D vs. l/rj s space in which the circles indicate those values 
of D and l/r/ s for which nonlinear saturation was achieved. The evolution is governed 
by the 3d electrostatic set of equations (0 EJ)- The other parameters have the same 
values as that chosen in Fig. 5. 

Fig. 7 Plot of the ratio r yrj s /k'^ vs. the permissible range of k z in simulation for rj s = 1.26 xlO -6 
(dotted line) and r) s = 1CT 5 (+ sign). 

Fig. 8 The slice plots for density and vorticity V 2 in the nonlinear state for two different 
values of rj s (viz. 10 -5 and 1.26 x 10~ 6 ). The other parameters and the governing 
evolution equation are same as that of Fig.5. 

Fig. 9 Scatter plots between n and tp and between V 2 <p and tp initially at t — and at the 
stage where saturation is achieved at t = 150. The data for n and (p is the same as 
that of Fig.5 with r) s = 10~ 5 . Clearly, theplots show no development of any kind of 
functional relationship unlike the 2d case [2|. 

Fig. 10 The numerical evolution of energy, streamer and zonal intensity per unit volume with 
time for the 3d electrostatic model of Eq. (|5l6|) ; for the Rayleigh Taylor (subplots in left 
column) and the resistive drift mode (subplots in the right column) for comparison. 
The plots in solid lines show the intensity in k z = modes and the dashed lines 
indicate the power in finite k z modes. 

Fig. 11 Slice plots showing the three dimensional density and potential structures in the linear 
(t = 150) and the nonlinear (t = 450) regimes for the resistively destabilized drift 
waves. 

Fig. 12 The evolution of total energy/volume for the case when L x /L y = 1 (solid line) and 
when L x /L y = 4. 

Fig. 13 The upper subplot shows the evolution of kinetic (solid line), pressure (dashed line) and 
the magnetic (dotted line) energy for the fully generalized electromagnetic simulations 
of Eqs.(0-HH). The lower subplot shows evolution of power in zonal (solid line) and 
streamer modes (dashed line). 
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Fig. 14 The slice plots showing three dimensional patterns for n, (J) and ip m the linear t = 30 
and the nonlinear t = 150 regimes for the electromagnetic simulation corresponding 
to the saturated state (r/ s V^| = 1) of Fig. 13. 
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